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Abstract We consider the constrained optimization of excitatory synaptic input pat- 
terns to maximize spike generation in leaky integrate-and-fire (LIF) and theta model 
neurons. In the case of discrete input kicks with a fixed total magnitude, optimal input 
timings and strengths are identified for each model using phase plane arguments. In 
both cases, optimal features relate to finding an input level at which the drop in input 
between successive spikes is minimized. A bounded minimizing level always exists 
in the theta model and may or may not exist in the LIF model, depending on pa- 
rameter tuning. We also provide analytical formulas to estimate the number of spikes 
resulting from a given input train. In a second case of continuous inputs of fixed total 
magnitude, we analyze the tuning of an input shape parameter to maximize the num- 
ber of spikes occurring in a fixed time interval. Results are obtained using numerical 
solution of a variational boundary value problem that we derive, as well as analysis, 
for the theta model and using a combination of simulation and analysis for the LIF 
model. In particular, consistent with the discrete case, the number of spikes in the 
theta model rises and then falls again as the input becomes more tightly peaked. Un- 
der a similar variation in the LIF case, we numerically show that the number of spikes 
increases monotonically up to some bound and we analytically constrain the times at 
which spikes can occur and estimate the bound on the number of spikes fired. 
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1 Introduction 

A major component of theoretical neuroscience is the study of how various neuronal 
models respond to synaptic inputs. Indeed, chemical synaptic transmission offers a 
specific mechanism for the encoding of information that an organism senses from 
the external environment, filtered by the internal state of the organism. The functions 
performed by particular neurons and neuronal networks are in part determined by 
the nature of the inputs that they receive and are in part a result of the responses 
they generate to these inputs, due to their intrinsic properties. Thus, understanding 
neuronal input-output transformations represents a centrally important scientific goal. 

Although the framework for incorporating synaptic inputs into computational 
models is well established, and the computational implications of such inputs have re- 
ceived significant attention, optimization problems involving synaptic inputs are not 
well represented in the literature. In this paper, we consider such a problem, namely 
what is the optimal way to tailor synaptic inputs, subject to certain constraints, to 
maximize the number of spikes that a neuron will fire? 

In fact, we consider two variations on this problem, one based on maximizing 
the total number of spikes fired and one focused on maximal firing within a pre- 
scribed time interval. There are several reasons that maximizing numbers of spikes 
may be a biologically relevant neuronal goal. Since neurons operate under conditions 
in which efficient resource use could be evolutionarily advantageous, it could be use- 
ful if, subject to some constraint on the amount of input that is available, the synaptic 
input time course could be tailored to achieve the largest possible number of spikes. 
Certainly, there are brain areas, including areas of visual cortex and somatosensory 
cortex, where it appears that intensity of firing encodes stimulus information, with 
neurons showing maximal firing under optimally preferred stimulus conditions [1-3]. 
Similarly, a sufficiently high firing rate within a given time window may be needed to 
overcome inhibition or to outcompete activity of other neurons to influence a down- 
stream readout neuron (for example, [4-6] and many more recent works). Which 
input time courses yield maximal spiking will depend on the intrinsic properties of a 
neuron or neuronal model, and characterizing optimal input time courses for models 
may provide information about which neuronal coding functions they are well suited 
to represent. 

In an earlier paper, a calculus of variations approach was used to determine the 
input current of minimal amplitude that causes a neuron to spike at a specified target 
time [7]. That work considered a phase model for a spiking neuron, with the evolution 
of phase x g [0, 27r) given by 

Jjc/Jr = /(jc) + Z(jc)/(0, 

where the impact of the input current I(t) is modulated by the phase sensitivity func- 
tion Z(x). Similarly, earlier work analyzed optimal weak inputs to start or stop repet- 
itive spiking in a general biological oscillator, with dynamics expressed in polar co- 
ordinates as 

dr/dt = €/(Osin(0), 

(1) 

d(l)/dt = 1 +€/(Ocos(0)/r 
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with € small [8]. Particular biological systems with dynamics qualitatively equivalent 
to those of (1) were also considered as examples to illustrate a calculus of varia- 
tions approach to this problem, involving numerical solution of a system of ordinary 
differential equations obtained through the introduction of Lagrange multipliers. To 
our knowledge, ours is the first work to focus on the optimization of inputs, subject 
to particular constraints, to maximize the number or rate of neuronal output spikes 
fired. 

To initiate research in this direction, we consider input optimization in two well- 
known mathematical models for single neurons, the leaky integrate-and-fire (LIF) 
model and the theta model [9-1 1], the mathematical formulations of which are given 
in Section 2. Both of these are scalar models, which allows us to use certain analytical 
approaches, rather than relying entirely on numerics, and we tune both models to 
be silent in the absence of inputs. One approach that we follow is to treat synaptic 
inputs as events with discrete onset times, which yields a phase plane representation 
of the co-evolution of an intrinsic neuronal variable and a synaptic input variable, as 
introduced previously for the LIF and theta models [12]. This approach would apply 
equally well in parameter regimes for which the models are intrinsically oscillatory, 
rather than silent, but we stick with the richer silent case. The LIF model was also 
considered by Borgers and Kopell, who proved (1) if a collection of phase- shifted, 
identical time-periodic inputs cause a spike in an LIF neuron, then the same inputs 
will also cause a spike if they are synchronized, and (2) if a constant input of size 
A causes a spike, then any time-periodic input with time average A will also cause 
a spike [13]. Both of these results show that the power of inputs to induce spikes in 
the LIF model increases as the input time course is made less uniform. Our results 
represent an extension of this idea, providing information about specific time courses 
that are optimal, not just for the generation of a single spike but for maximizing spike 
outputs. While the LIF model is a reasonable representation for a passive neuron, the 
theta model can be rigorously derived as a normal form for Type I spiking neurons [9- 
11], which feature a transition from silence to oscillations through a SNIC bifurcation 
[14]. Thus, consideration of the theta model allows us to explore how our results 
extend to a case with additional biological relevance. Interestingly, although both 
models can fire spikes at arbitrarily low frequencies, differences in their responses to 
inputs have been noted in previous work [12], which further motivates us to continue 
to compare them here. 

In addition to discrete inputs, we also consider a continuous input formulated so 
that a particular parameter shapes the input time course and ask how that parameter 
should be tuned to achieve the maximal number of spikes within a given time win- 
dow. The absence of a threshold and reset in the theta model is convenient for the 
application of variational techniques in this case. For the LIF model, we do not ap- 
ply such techniques, but we nonetheless derive some results about the influence of 
the input shape parameter on the spike train that results from the input. We shall see 
that certain properties of each model observed in the discrete case carry over to the 
continuous case, yet there are some differences worth noting as well. 

The remainder of the paper is organized as follows. In Section 2, we analyze the 
LIF model with discrete synaptic kicks, the cumulative sizes of which are constrained. 
We introduce phase plane structures that are useful for analysis of the resulting op- 
timization problem and discuss some example strategies for controlling the timing 
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of inputs before moving on to prove some results about these structures and optimal 
strategies. In particular, we show that two different scenarios are possible, depending 
on model parameters, and these yield different optimal input strategies. We also show 
how analytical estimates for the number of spikes fired can be derived for any given 
input time course. In Section 3, a somewhat similar analysis strategy is applied to the 
theta model with discrete synaptic kicks. We find that the phase plane structures for 
the theta model resemble one of the possible cases for the LIF model, with a cor- 
responding similarity in optimal strategies. In Section 4, we turn to the analysis of 
continuous input for the theta model, seeking a maximal number of spikes in a fixed 
time window. We define the continuous input so that its shape can be modulated 
by a parameter yet its integral over all positive time is independent of the value of 
that parameter. Variational techniques yield a boundary value problem that we solve 
numerically to find locally and globally optimal parameter values. In Section 5, we 
consider the continuous input optimization problem for the LIF model, for which 
the discontinuous reset prevents application of the same variational techniques. We 
provide direct analysis showing that all spikes must be fired within a particular time 
interval and characterize the behavior of this interval as the input parameter is varied. 
Moreover, we prove that the number of spikes saturates as this parameter increases. 
Some similar analysis for the theta model, showing that spiking is lost as this param- 
eter increases, appears in the Appendix. We conclude with a discussion in Section 6, 
where we summarize our results and discuss the novelty and relevance of our findings 
in the neuroscience setting. 



2 Integrate-and-fire (LIF) model with synaptic kicks 

Model 

We consider the dynamics of an LIF model neuron with excitatory synaptic inputs as 
governed by the equations 

v' = 1 -v-g{v-E), (2) 
g' = -fig, (3) 
together with the reset condition 

v(r) = Vth =^ v(t^) = Vr. (4) 

Equation (2) can be derived from a conductance-based equation Cv^ = I input — 
J2j gj(^ ~ Ej) ~ gsyn(V — Esyn) with fixed intrinsic current conductances gj , but we 
think of it as a nondimensionalized abstract model in which voltage intrinsically con- 
verges to a baseline / and E is the reversal potential of a synaptic input with strength 
^ > 0. We assume Vr < I < Vth, such that no spikes are fired in the absence of in- 
put, and E > Vth, and we consider the invariant half-plane {g > 0} within the (i;, g) 
phase plane, where (/, 0) is the unique stable critical point. Further, we represent the 
excitatory input by the equations 

g(tn)=gOn)+kn, (5) 
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(6) 



n=l 



foxkn G (0, G], ^ = 1, 2, . . . , with N >\ finite and G > 0 fixed in R. Equations (3) 
and (5) show that each input kick can be chosen to arrive at any time and instanta- 
neously updates the value of g when it arrives and that the synaptic conductance g 
always decays exponentially between kicks. Equation (6) states that the sum of all 
inputs, however they may be divided, is always equal to a fixed input allowance G. 
In the subsequent subsections, we assume that I , E, Vth, Vr, and are fixed and con- 
sider how to partition and time the input G to yield the greatest number of threshold 
crossings, or spikes. Without loss of generality, we take Vr = 0. First, we discuss a 
phase plane representation of this problem and consider some example strategies. 

Phase plane structures and basic strategies 

We illustrate some key structures in the phase plane for system (2), (3) in Figure 1. 
The i;-nullcline, based on equation (2), is the curve {g = (v — I) / (E — v)} . Denote the 
trajectory through the point := (vth, g^ = (vth — n/(E — Vth)) on this curve by 
Fo, which is tangent to the line {i; = Vth}- Fq partitions the set {(v, g) \ v < Vth, g > 
0} into a set of initial conditions that yield spikes, namely those above Fq, and a 
set of initial conditions that do not yield spikes, namely those below Fq. Let F^ 
denote the intersection of Fq with {i; = i;^} and let g^ denote its g coordinate. The 
point F^ := (vth, g^) is on the threshold line, and there is a trajectory Fi that flows 
through this point. The band between Fq and Fi, with Vr < v < Vth, consists of the 
set of initial conditions from which trajectories yield exactly one spike. Denote this 
band by Bi . Similarly, for each natural number n, we can define a curve F^, such that 
each trajectory with an initial condition in the band Bn between F^_i and F^ yields n 
spikes. Note that this band structure exists if parameters are altered such that I > Vth, 
although the minimal band is shifted to negative g values, and hence the methods that 
we discuss easily generalize to this case. 

Fig. 1 Key structures in the J? a , ^ 



phase plane for system (2), (3), 
Note that G, Vfh, and / are 
given parameters, while the 
other structures depend on the 
system dynamics. 




V 
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The band structure that characterizes sets of initial conditions that yield particular 
numbers of spikes calls to mind several natural strategies for doling out input kicks 
to maximize spike output: 

Critical kicks 

Once a trajectory is below Fq, it approaches (/, 0) asymptotically. Let Yc denote the 
^-coordinate of the intersection {v = /} H Fq. One possible strategy is to give an 
initial input of size yc as well as subsequent inputs of size yc each time the trajectory 
reaches a sufficiently small neighborhood of (/, 0); see Figure 1. This critical kicks 
strategy yields n spikes where nyc < G < (n -\- l)yc. 

Big kick 

A possible disadvantage of the critical kicks strategy is that the increment yc to 
achieve a spike may exceed the width y^ between bands F^ and F^_i for ^ > 1. 
To ensure that this increment is only encountered once, taking inspiration from the 
power of synchronized inputs [12, 13], another reasonable strategy is to give a single 
big kick of size G, all at once, to the resting cell (Figure 1). 

Reset and kick 

An input of size yc is sufficient to push the voltage across the threshold for single 
spike initiation. In the big kick strategy, the additional available input G — is pro- 
vided together with the yc. It is possible that this additional input could deliver more 
spikes if it were delivered separately from the initial yc. We can define a strategy, 
for example, in which an initial kick of size yc is given to elicit a spike. As soon as 
this spike is fired, the cell is reset and the remaining input allowance G — is given. 
Clearly, this reset and kick strategy would make sense if the bands Bn were narrowest 
at i; = 0. 

Threshold kick 

At the other extreme, the bands Bn might be narrowest at i; = Vfh- In this case, a 
possible optimal strategy would be the threshold kick strategy, defined by giving the 
initial kick of size yc and following this with a kick of size G — yc just before thresh- 
old crossing occurs (Figure 1). 

Intuitively, it is reasonable to think that if ^ is large, such that inputs rapidly decay, 
then it makes sense to dole out inputs in minimal pieces, such that something like the 
critical kicks strategy may be optimal. Alternatively, if ^ is small, such that inputs 
decay slowly, then it makes sense to make inputs available as early as possible, such 
that one of other strategies is likely to be optimal. To analyze more carefully which 
strategy is optimal, it will be helpful to define a band width 5^ (f) as the distance from 
Tn-i to F^ in the ^-direction for each fixed v G [0, Vth\- With this definition in hand, 
we note that 5^+1 (i;^/^) = 5^ (i;^) for /t > 1, and hence the reset and kick and threshold 
kick strategies are effectively the same strategy, yielding the same number of spikes 
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(Figure 1). We also let 8oo(v) = lim^^^oo ^niv), v ^ [0, Vth], if this limit exists. Very 
roughly speaking, the critical kicks strategy will yield approximately G/yc spikes 
while the other strategies will induce about G/8oo(v) spikes for some v, at least if 
8n(v) converges to 8oo(v) quickly. Thus, comparison of yc and 8oo can be used to 
give an initial suggestion of what strategy to follow. 

The value of yc can be observed numerically by backwards integration from 
('^th, Sq) until V = I. Alternatively, it may be optimal to replace yc by the distance 
j/^ := ^~ — as can be computed by backwards integration from (vth, Sq) up to 
(Vr, Sq), and give kicks of size yc after each spike, after the initial kick of size y^; 
we will still refer to this as a critical kicks strategy. 

An approximate value of 8oo(v) can be derived as follows. From (2), (3), the slope 
s of the vector field at any point in the phase plane is given by 



For V G [Vr, Vth], 



s(v,g) = j ^7 -• (7) 

I -v-g(v-E) 



^ 

s'^iv) := lim s(v, g) = < 0. 

^^oo E — V 



The magnitude of the change in g over one spike cycle is 



rvr 

:= / s^dv 
= / -^r^. dv = l3ln 

J V, 



Vth - E 



(8) 



By construction, 8oo(Vr) = 8oo as given by (8). But since the value of g at reset for 
the trajectory forming the upper bound of one band is the value of g at threshold for 
the trajectory forming its lower bound, and we have taken the limit as n ^ oo, we 
can also estimate 8oo(vth) using (8) and indeed, using similar translation arguments, 
we estimate 8oo(v) = 8oo for Sillve[vr, Vth]- 

Comparison of yc (or yc) and 8^ suggests whether or not the critical kicks strategy 
will elicit more spikes than the other strategies we have described. If not, then we need 
additional arguments to assess the relative effectiveness of these alternative strategies. 
In fact, regarding alternative strategies, we have the following result: 

Proposition 1 The big kick strategy always yields at least as many spikes as the reset 
and kick {and equivalently, the threshold kick) strategy. 

Proof The reset and kick strategy yields m + 1 spikes, where m is the largest integer 
such that 



^8n(vth) <G-yc. 



n=l 
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Using equation (7), compute 

ds - v) 



dg (I-v-g(v-E))^ 



(9) 



We can see from equation (9) that if v < I , then the slope s becomes more negative as 
g is increased, and if i; > /, then the slope s becomes less negative as g is increased. 
Thus, the bands are narrowest at v = I; that is, 8n(I) < mm{8n(vr), 8n(vth) = 
8n-\-i(Vr)}. Hence, the big kick strategy, which elicits + 1 spikes for the largest 
integer such that 

rrib 

J28n(I)<G-yc, 

n=l 

always generates at least as many spikes as reset and kick. □ 
Band width estimation 

In the previous subsection, we introduced a small number of intuitively reasonable 
strategies for eliciting the maximum number of spikes from model (2), (3) using 
a constrained input. We also used a phase plane approach to define a natural band 
structure, along with a corresponding idea of a band width 8n(v), which we used to 
show that two of these, the reset and kick and threshold kick strategies, will never be 
optimal. This structure can also be used to obtain an intuitive idea of which conditions 
favor a big kick strategy and which conditions favor the critical kicks strategy of giv- 
ing many small kicks of the same particular size. Next, we use some approximations 
to derive additional quantitative information about 8n that can be used to determine 
more globally the optimal input strategy. Henceforth, in addition to assuming that 
= 0, we for convenience set vth = with 0 < / < 1 < 

We can estimate the magnitude 8(g) of the change in g that occurs over one spike 
cycle using the slope s(v, g) given in equation (7), 



/ I-v 



^(8)= I -r ^7 Trdv 

g(v — E) 

^ (10) 
_ fig I + gE-l-g 

^ m , 

1+g I+gE 

where we have approximated g by a constant to estimate the integral. Note that 
for each n, the band widths 5^(0) = 5^+1 (1) from the previous subsection are ap- 
proximately equal to 8(g) for certain corresponding choices of g; for example, 
8i(0) ^ 8((gQ + g^)/2), where Fi intersects {v = 0} at (0,g^). More generally, 
it is not necessary to choose a g associated with the boundary of a band, as defined 
from the previous subsection, in order to compute 8 (g) . 

We can investigate the spikes of a cell by analyzing (10). This approach yields the 
following result. 
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Proposition 2 If E -\- 1 — 2EI >0, then 8(g) is a monotone decreasing function ofg. 
If E -\- 1 — 2EI < 0, then 8(g) has a unique local minimum at a finite, positive value 
g = go- 
Proof Calculating the derivative of 8(g) with respect to ^, we have 

d8(g) 



where 

/(g):=ln 
Furthermore, 



dg (1+^)2"^^^^' 
I + gE \ g(E-I)(l + g) 



I^gE-l-gJ (I^gE-l-g)(I^gE) 

df(g) ^ (E- I)(l + g)(g(E + I - 2EI) + 2/(1 - /)) 
dg (i^gE-l-g)Hl^gE)^ 



(11) 



Equation (11) shows that E -\- 1 — 2EI is indeed a key quantity. 

Suppose now that E + I - 2EI > 0. If £ + / - 2EI > 0, then df(g)/dg > 0, 
such that f(g) increases as g increases. Define 

E E-I 

fi = liminf f(g) = In . 

g^+oo-^^^' E-1 E(E-l) 

Since 

dfi E + I - 2EI 

je = -eW^-'^ '''' 

fi increases as E increases. As 

liminf /i = 0, 

we have f\ < 0. Therefore, f(g) < 0 for all g and hence d8(g)/dg < 0. In another 
words, under the original approximation used to obtain (10), 8(g) decreases, and thus 
there is less change in g across each cycle from reset to threshold, as g increases. 

Next, suppose that E -\- I — 2EI < 0. Under this condition, df/dg changes 
signs, with df/dg > 0 for ^ g (^, 2ei-~e-i ^ df(g)/dg < 0 for ^ g 
^ 2Ei^E-i ' +^)- From expression (12), we have dfi /dE <Q and liminf^^+oo fl = 
0, such that /i > 0 for all E and /(^) > 0 for g e ( 2ei-~e1i ^ +^)- Thus, there ex- 
ists a unique point go such that f(g) < 0 for all g e (-^E^, go) and f(g) > 0 for all 
g e (go, +oo). Correspondingly, 8(g) decreases when g e (-^E^, go) and increases 
when g e (go, +oo), where go is the zero point of f(g) = 0, and the proof is com- 
plete. □ 

In the case of £" + / — 2£'/ > 0, the monotonicity of 8(g) suggests that for a 
trajectory evolving from an initial condition of (v, g) = (0, ^(0)) to a final condi- 
tion (v, g) = (I, g(t)), the drop ^(0) — g(t) should be smaller for larger ^(0). From 
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Fig. 2 The red lines display the estimate of 8{g) for (2), (3) obtained by our calculations, while the blue 
lines show 8(g) from direct simulation for E = 1.2, I = 0.7, and (a) ^ = 0.5, G = 10 or (b) /3 = 2, 
G = 20. 



Figure 2, we can see that, while the approximation used to derive equation (10) intro- 
duces an error in relative to the actual change in g computed from direct simulation 
of trajectories with E -\- 1 — 2EI >0, the error appears to be small and the mono- 
tonicity of 8(g) appears to be correct. Similarly, a numerically computed example of 
8(g) for parameters that yield £" + / — 2£'/<0is shown in Figure 3. 

Optimal strategies and spike counts 

Based on the previous two subsections, we conclude that if E -\- 1 — 2EI >0, then 
the loss of input with each spike, measured by 8(g), decreases as g increases and, 
furthermore, an input of a fixed size will cross the most spiking bands if it is given at 
V = I. Hence, of all possible strategies for eliciting spikes with an input of total size 
G, the one that yields the most spikes is what we earlier called the big kick strategy, 
unless Yc (or Yc) is sufficiently small that avoiding the bands altogether by following 
the critical kicks strategy is optimal. The optimal strategy when £" + / — 2£'/<0is 
to provide a kick that puts g at approximately the g value where the minimum of 
8(g) occurs, and then provide as many kicks as possible of size 5 (^o). again assuming 

Fig. 3 A numerically 
calculated example of 8(g), with 
a unique local minimum, for 
parameter values E = 2, 
I = 0.7, ^ = 0.05, G = 10, such 
that £ + / -2^7 <0. 
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that Yc, Yc are above a certain size. We will next perform some additional calculations 
that can provide estimates of numbers of spikes resulting from any strategy, which 
can be used with a minimum of calculation to compare the results of particular input 
sequences. 

We first suppose that E-\-I — 2EI <0 and consider a generalization of the optimal 
strategy described above. That is, we assume that an initial input Gi > go is given and 
then, once g evolves to some neighborhood of kicks of size 8 (go) are repeatedly 
applied until the remaining input falls below 8 (go). We now estimate the numbers of 
spikes ^2 fired for each strategy of this type. Let Qi denote the number of spikes fired 
during the initial time period when g drops toward ^o, let denote the number of 
spikes fired during the final time period after the available input is depleted, and let 
denote the number of spikes fired during the intervening period when repeated 
kicks of size 8 (go) are given. Clearly, 

G - Gi 

Q3 ^ (13) 

8(go) 

and ^2 = ^2i + + ^^3, so it remains to estimate ^2i and ^22. 

Because of the shape of the function 8(g), the largest 8(g) during the initial time 
period will be associated with the first spike fired, while the largest during the final 
period will be associated with the last spike fired. We can estimate the drop 8(Gi) 
in g up to the firing of the first spike from equation (10). To make this estimate 
relevant to other spikes early in the spike train, we take i;(0) = 0 rather than i;(0) = / . 
Approximating the level of g in equation (10) by := G/ — 5(G/)/2, we obtain 

S(G^)^Sr.= -^ln '-';f-'^'\ (14) 

l-\-gi / + Egi 

Next, we estimate 8(g) for the final spike fired, which we call 8f. To do this, we 
assume that when the final spike is fired, the trajectory reaches the lower bound on 
the g values that can yield a spike, namely the point of intersection of the u-nullcline 
and {i; = 1}, at which ^ = ^^ = (1 — /)/(£" — 1). We also use the intermediate value 
of g across the trajectory ro(i'), namely (I — I)/(E — I) -\- 8f/2, as the value of g 
for equation (10), which yields 

^ P[8f/2 + (l-I)/(E-l)] 

^ l+5//2 + (l-/)/(£-l) 

^^^ I-l + (E- l)[8 f/2 + (1 - I)/(E - 1)] ^^^^ 
I^E[8f/2^(l-I)/(E-l)] 

Now, to obtain an estimated spike count as g decays from G/ to we approxi- 
mate 8(g) over each spike by the average of its two extreme values, (8 (go) + 8i)/2. 
This approximation yields 

^ ^ G,-(go-.(go)/2) 
iS{go)+Si)/2 

^ Springer 



Page 12 of 35 



Wang et al. 



Similarly, once the input is used up, spikes continue to be fired as g decays from 
to approximately (1 — /)/(£ — 1), and the number of additional spikes that result is 
estimated by 

^ ^ go + 8(go)/2-(l-I)/(E-l) 
(8(go)^8f)/2 

In the above equations, we have taken into account that the trajectory may not be 
reset precisely at go but rather somewhere within an interval approximated by (^o — 
8(go)/2, go + 8(go)/2). Because this may lead to an overestimation by one or two 
spikes, we set + = ^1 + — 1. The total number of spikes fired is finally 
estimated by 

^2 = ^21 + ^22 + ^^3 (18) 

using equations (13)-(17). 

The calculation can be easily generalized for input patterns that push g above and 
below go multiple times, although they will be non-optimal by our earlier arguments. 
Similarly, for an initial kick < go and g < go for all time, the smallest 8(g) avail- 
able is 8(Gi) and we can estimate the number of spikes resulting from partitioning 
the input into kicks of size 5(G/) by the equation 

G-Gi Gi -(1 -/)/(£-!) 

^2= ^ + ^ ^ ^-^ -. (19) 

8i {8i^8f)/2 

If E -\- 1 — 2EI >0, then the same calculations still apply. The big kick strategy 
is optimal here, yielding a number of spikes estimated by 

(8(G)+8f)/2 

Generalizing, a strategy of giving an initial input G/, followed by repeated kicks of 
size 8i until the input is depleted, yields a number of spikes estimated by 

G - Gi Gi - (1 - I)/(E - 1) 

^2= - + — ^ — (21) 

8i (8i+8f)/2 

Figure 4 shows comparisons of our spike estimates from equations (18) and (21) with 
numerical computed counts of spikes, illustrating that our estimates can be reason- 
able. 

In fact, in the case of £ + / — 2£'/ > 0, we underestimate the number of spikes 
fired for large G and G/. This underestimation results because we average 8(G) or 
8i with 8f in the denominator of equation (20) or (21), whereas most spikes yield 
decreases in g that are much smaller than 8 f . Improved estimates in such cases be 
obtained by weighting this denominator more toward 8(G) or 8i, which will decrease 
the denominator and thus will always yield predictions of additional spikes for larger 
kicks, relative to the formulas in equations (20), (21). An example resulting from 
extreme weighting, replacing the average of 8i and 8 / with 8i alone, is also shown in 
Figure 4, as is a similar example for the case of E -\- 1 — 2EI < 0. 
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(a) (b) 

Fig. 4 Number of spikes fired as a function of G/ for G = 100. The blue points show results from di- 
rect numerical simulations of (2), (3) while the red and black points show the predictions from equa- 
tions. Specifically, the red points are derived from (a) equation (21) with E = 1.2, I = 0.7 such that 
E + I - 2EI > 0 and (b) equation (18) with E = 2,I = 0.1 such that E + I - 2EI < 0. The same param- 
eter values were used for the black points but the average iS(gi) + 8f)/2is replaced by S(gi) in equation 
(21) in (a) and the average {8{gQ) + 8{gi))/2 is replaced by 8{gQ) in equation (16) in (b). 



In summary, equations of the form (18)-(21), each requiring calculation of only 
a small number of quantities, can be used on a case by case basis to estimate the 
numbers of spikes that will result from a given strategy and therefore to compare 
strategies. These formulas provide for an informed comparison between the two types 
of big kick strategies determined to be optimal for the two distinct cases of £" + 
I — 2EI > 0 and E -\- I — 2EI < 0, respectively, and the critical kicks strategy. 
Furthermore, now that we have defined 8(g), we can give a more precise variation 
on the calculation of equation (9) made in the subsection on phase plane structures 
and basic strategies to show that truly optimal strategies (other than the critical kick 
strategy based on Yc) provide kicks with v = I, so each strategy should include a 
time shift so that kicks are given when this condition is met, rather than with i; = 0 
or i; = 1. Specifically, if > g2, then for i;o G [0, 1], 



^O) = / 



8(gi,vo)-8(g2,vo)= I —dv 

gi(v-E) 



i 



dv. 



Jo I -V- g2(v - E) 

Calculating the derivative of the above equation with respect to vo yields 
d{8{gi, vo) - 8{g2, vo)) P(I - vo)(gi - g2) 



dvo (I -vo- glVo + glE)(I -vo- g2V0 + g2E) 

a quantity that is positive for vo < I and negative for vo > I . Hence, the additional 
input needed to cross bands is minimal for kicks given aX vo = I , agreement with 
Proposition 1 . Finally, it is not difficult to see from examination of the above spike 
counts and equation (10) that, with other parameters fixed, increases in yield fewer 
spikes, as expected from the corresponding faster decay of g, while increases in E 
and / yield more spikes, as expected from the increased rate of change of v. 

^ springer 



Page 14 of 35 



Wang et al. 



3 Theta model with synaptic kicks 



Model 



We next consider a theta model neuron receiving positive synaptic excitatory kicks, 
governed by the equations 



0' = l-cosO-\-{b-\-g)(l+ cosO), 




where Z? < 0 is a parameter. We consider 0 e [— tt, tt] mod 27r, and the neuron is 
said to fire when 0 increases through tt and is effectively reset to —jr. With ^ = 0, 
corresponding to the absence of excitatory inputs, and Z? < 0, which is the case we 
consider, the theta model (22) has two critical points, namely a stable fixed point at 
Os = — arccos < 0 and an unstable fixed point at Ou = arccos > 0. Moreover, 
we assume that the arrival of, and constraints on, synaptic excitation are identical to 
those for the LIF model, given by (5), (6), with G > 0 fixed. Note that, as in the LIF 
case, everything that we do in this section would still apply if we chose b sufficiently 
large that the model fired spikes in the absence of input, but we stick with the Z? < 0 
case to include the additional effects associated with the requirement of a minimal 
input for spiking to occur. 



Existence of an optimal g 



We proceed by approximating the amount by which g decreases as a trajectory 
evolves from (— tt, gi) to (tt, g /), analogously to our calculations in Section 2. Fixing 
g at some value between gi and ^/ , we have 



'">=r i- — 



cos6> + (/? + ^)(l + cos6>) (23) 
pgjc/^bTg. 



A straightforward calculation yields the following result. 

Proposition 3 The function 8(g) defined in (23) has a unique local minimum at g = 
-2b. 



Proof We calculate 

d8(g) P 



zir — 



dg 2(/7 + ^)2/3^ 

fi(2b^g) 



2(b + g)y- 



:7t. 



So d8(g)/dg < 0 when g < —2b, and d8(g)/dg > 0 when g > —2b. Clearly, 8{g) 
has the minimum value 2pTt^J^-b at g = —2b. □ 
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Fig. 5 The red lines represent S{g) for the theta model from equation (23), while the blue lines give the 
size of the change in g observed from direct simulation of system (22) for (a) b = —2, ^ = 0.05, G = 10 
and (b) b = -20, ^ = 0.3, G = 100. 



In Figure 5, we validate the approximation of holding g constant in equation (23) 
by showing that such a minimum exists in direct simulations of the theta model. 

Optimal strategies and spike counts 

We will now estimate the number of spikes that will result from a given input alloca- 
tion strategy. To make effective estimates, it is helpful to estimate the minimal value, 
call it g, such that a trajectory starting from (—7t,g) will result in a spike if and only 
if ^ > ^. To do this analytically, we seek g such that the trajectory from (—7t,g) 
reaches (0, —b) and thus crosses the ^-nuUcline and converges to Os', see Figure 6. 
Although there are other trajectories with initial g values above this g that also con- 
verge to Os, instead of spiking, by crossing the ^-nullcline at points with ^ > 0 and 
g < —b, this approach nonetheless gives a reasonable approximation to the true value 
ofg. 

We again approximate g as a constant over the trajectory of interest, namely g = 
(g — b)/2, which is the average of the initial value of g and the value of g at ^ = 0 



Fig. 6 Schematic illustration of 
the (6, g) phase plane including 
our approximation to g and g. 
The shaded region indicates the 
set where 6^ < 0. 82 is the 
change in ^ as ^ evolves from 
—TT to 71 for the firing of the 
neuron's final spike. 




last spike 
fired 



-K (65,0) 9=0 (eu,0) n 
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(Figure 6). This approximation yields 

g+b= r — do, 

^ Jj, l-cosO + (b + g)(l+cosO) 

and thus 

g + b = — , ^ 

V(| + Z7)/2 

which provides an impHcit estimate for g. 

Now, if g > —2b, then all inputs that yield spikes must push the trajectory to g 
values for which d8(g)/dg > 0 holds. Thus, pushing g to progressively larger values 
yields fewer spikes, and a critical kicks strategy, with initial kick size g and subse- 
quent kick sizes 8(g), is optimal, yielding approximately 

G-g 

Q = — ^ (24) 
8(g) 

spikes. 

If ^ < —2b, then there is a tradeoff: investing more input initially, up to about size 
—2b, will push the trajectory to a region where 8(g) is minimal. On the other hand, if 
less input is initially invested, then there is more input remaining to give subsequent 
kicks. We ignore strategies in which an initial input Gi > ^ is given, some nonzero 
number of spikes is fired, and then an additional large input spanning multiple spiking 
bands is given, since these can be seen always to be non-optimal. Suppose first that 
the initial input has a size Gi < —2b. Given the shape of 8(g), the optimal strategy is 
to expend the remaining input on 

G-Gi 

(25) 

8(Gi) 

kicks of size 8(Gi), analogously to the strategy behind equation (24). Which Gi < 
—2b is optimal depends on the sizes of g,b, and 8(g). 

Alternatively, the other possible optimal strategy, if G > —2/?, is to take Gi >—2b 
and to try to maintain g values close to g = —2b. For an initial input Gi > —2b, we 
estimate the number of spikes for such a strategy using a similar approach to that used 
in Section 2.4, breaking up the estimate into an initial period of decay of a period 
of kicks to keep g near —2b, and a final period of spiking until g drops below g. Once 
the initial input is given, 8(g) for the first spike is approximated by a solution to the 
equation 

^(Gi-8i/2) 
8i = 7t. 
^b^Gi -8x12 

We will also use an estimate of the 8(g) for the last spike fired (Figure 6), obtained 
from the equation 

)g(^2/2 + g) 
82 = =7t. 

^b^82/2 + g 
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Now, for g < —2b, recall that d8(g)/dg < 0 when g < —2b and d8(g)/dg > 0 
when g > —2b. Thus, during the period when g > —2b, the largest 8(g) is given by 
8i, and while g < —2b, the largest 8(g) is 82- The smallest 8(g) is of course 8(—2b). 
Following our earlier strategy of estimating 8(g) by the average of its largest and 
smallest values, and approximating g over one spike interval by the initial value mi- 
nus half of the drop in g that occurs during that interval, we obtain our estimated spike 
counts. Specifically, during the initial period of input decay from G/ to approximately 
—2b, our estimate is 

^ _ Gi-8(-2b)/2^2b 
^~ (8(-2b) + 8i)/2 

During the final period of input decay from approximately —2b to approximately g, 
our estimate is 

^ ^ _2b^8(-2b)l2-g 
^ (8(-2b) + 82)l2 

Finally, the number of repeated spikes from the critical kicks of size 8(—2b) given 
until the remaining available input G — Gt is depleted is about 

8(-2b) 

Given that and could each overestimate the number of spikes by one, we 
estimate the total number of spikes with input G as 

^2 = ^2i +^22 + ^23 - 1. (26) 

A similar estimate can be obtained from other patterns of inputs. 

In summary, we have two candidate optimal strategies when g < —2b. Depend- 
ing on the relative sizes of g and —2b and the shape of 8(g), it may be optimal to 
choose initial input < —2b and give repeated kicks of size 8(Gi) or it may be 
optimal to choose G/ ^ —2b -\- 8(—2b)/2 and provide repeated kicks of size 8 (—2b), 
in both cases repeating the kicks until the input is depleted. Figure 7 shows examples 
illustrating that the optimal G/ is indeed very close to —2b. 

We can extend this analysis one step further by determining the best value of 0 at 
which to deliver the input kicks. 

Proposition 4 Given an initial condition (—7t,g) and an input of a fixed size, the 
maximal number of spikes subsequent to input delivery is attained when these inputs 
are given with 0 = arccos . 

Proof Suppose ^1 > g2- We have 
Kgi,Oo)-8(g2,OQ) 

Pgl Pg2 



do. 



-cosO - (b-\- gi)(l-\-cosO) 1 -cos^ - (b-\- g2)(l-\-cos0)_ 
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°o (a) % (b) 

Fig. 7 Spike counts for the theta model (22). The red curve denotes the number of spikes estimated from 
equations (25), (26), while the blue lines show numbers of spikes observed in direct simulations of the 
theta model with G = 100 and (a) b = -5, P = 0.2; (b) b = -20, y0 = 0.1. For the simulations, an initial 
kick of size G/ was given. If G/ > —2b, then the trajectory was allowed to decay until g ~ —2b — 8(—2b) 
and then the remaining input was expended on kicks of size S{—2b). If G/ < —2b, then after an initial 
spike was fired, the remaining input was expended on kicks of size 8{Gi). These same strategies were 
assumed for the analytical estimation. 

Calculating the derivative of the above equation with respect to ^o, we have 

d(8(guOo)-8(g2,Oo)) 
dOo 

^ P(l^b^(b-l)cos0o)(gi-g2) 

[1 - cos6>o - (/? + + cos6>o)][l - cos6>o - (/? + g2)(l + cos6>o)] ' 

As d(8(gi,Oo) - 8(g2, Oo))/dOo > 0 for ^ [-^r, - arccos ]^) U (arccos {3^, Jr] 
and d(8(gi, vq) — 8(g2, vo))/dvo < 0 for ^0 ^ (—arccos arccos {3^), the dif- 
ference 8(gi,Oo) — 8(g2,Oo) will be maximal at ^0 = — arccos and minimal at 
^0 = arccos . Hence, a maximal number of spikes is attained from an input given 
with 0 = arccos . □ 



4 Theta model with continuous input 

In the analysis we have done so far, the input to the neuron arrives as a series of 
discrete kicks. An excitatory postsynaptic potential evoked by an individual input 
may have a more gradual rise, however. We now switch gears and consider how such 
an individual input, arriving with a continuous time course, can be optimally tailored 
to evoke a maximal response. This analysis is also relevant to a situation in which 
a neuron receives input from a very large presynaptic population that fires in near 
synchrony, but with some spread in recruitment times. 

The theta model with a continuous input can be described by the equations 

= 1 - cos6> + (/? + y(t))(l + cos6>), (27) 
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with 

y(t) = Afi^tQxp-^\ (28) 

where b <0 and A, ^0 > 0 are parameters. The form of y(t) in equation (28) is often 
used computationally and has been specifically selected to ensure that its integral over 
the positive time axis is fixed at A for all values of fi. Unlike the case where the 0 
model received synaptic kicks, we now consider ^ g R, with a spike fired whenever 
0 increases through rut for an integer n. 

The question that we now address is, given a fixed set of intrinsic parameters b 
and A and fixed P > 0, how should p be selected to yield the maximum number of 
spikes within the time interval [0, P]l Note that when 0 = nn , 0^ > 0, so it suffices 
to find to maximize 0(P). 

Boundary value problem 

We will find optimal values of by numerically solving a boundary value problem. 
For numerical purposes, it is convenient to map t e [0, P] to 5- g [0, 1] using s = t/P. 
Correspondingly, let 'denote d/ds, such that equation (27) becomes 

^ = p(l - cos6> + (/? + y(s))(l + cos 6*)). (29) 

Next, differentiate 0 with respect to to obtain 

0^ = p{e^ sin^ - smO{b + y{s)) + y^{s){\ + cos^)), (30) 

where y^{s) = ^y{s) is given by 

y^{s) = ApPs exp-^^"(2 -pPs). (31) 

To the 0 and 0^ equations, we append the additional equations 

^=0, (32) 
/ = P. (33) 

Given the system of four equations (29), (30), (32), (33), we need a set of four 
boundary conditions. First, we set ^(0) = Os = — arccos((l +/?)/(! — /?)), so that 
the model neuron is at a stable rest state when input starts to arrive. Since this spec- 
ification is independent of yS, we have ^^(0) = 0. To find an extreme (with respect 
to fi) value of 0(P), we set 0^(P) = 0 as well. Finally, we take ^(0) = 0. We solve 
this boundary value problem (BVP) numerically, using XPPAUT [15], to obtain the 
optimum for any fixed P, A,b. 

Results 

If A and P are fixed, then varying yields different shaped input functions y(t), 
as illustrated in Figure 8. For fixed A, if y0 is sufficiently large, then the input is 
sufficiently concentrated in time that P becomes irrelevant for the number of spikes 
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Fig. 8 y{t) for A = l, P =2, 
md^ = l (solid), 2.316 
(dashed), and 4 (dotted). 
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fired. Alternatively, for smaller input arrives more gradually and P becomes rele- 
vant. These relationships are evident in Figure 9A, which shows a curve of solutions 
to the optimization BVP described in the previous subsection, plotted in ^ versus P 
space. Direct numerical simulations show that in fact the same 0(P) arises for other 
large as occurs for the value along the upper branch of the curve in Figure 9A, 
as long as P is sufficiently large that this curve is flat (Figure 10). Along the rest of 
the curve, the values shown truly represent local extrema. In particular, the lower 
branch of the curve represents local maxima for 0(P). For P values such that the 
upper branch is flat, the lower branch appears to actually consist of global maxima 
(Figure 10), while the upper branch represents local minima for sufficiently small P 
that this this branch has curvature (for example, P near 2), and for these P, there may 
be additional maxima for 0(P) at large fi, not shown in Figure 9A. 

For other values of A, the situation is similar to Figure 9 A; however, some sub- 
tleties do emerge. As shown in Figure 9B, additional local extrema of 0(P) may 
arise at small fi, yielding an interval of P values with three such extrema at small 

bracketed by two fold bifurcations. For example, for A = 8 and P = 10, there are 
a local maximum near = 0.31, a local minimum near = 0.57, and another local 




10 



15 

(B) 



Fig. 9 Optimal ^ versus P for (A) A = 7 and (B) A = 8. 
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Fig. 10 Time course of 6 for 
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P = 4, as plotted in Figure 9 A, 
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are ^ ^ 0.95 (lower branch) and 
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^ 7.28 (upper branch); 


e 3 
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identical (6 ^ 5.04) for all ^ 
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sufficiently large. 
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maximum near = 0.72, with corresponding 0(t) shown in Figure 11. Note that, as 
with A = 7 as well as other values of A, 0(P) saturates for sufficiently large fi. The 
folding structure in Figure 9B arises only for a small interval of A values that is diffi- 
cult to pin down precisely. It is also illustrative to view what happens to the families 
of B VP solutions as A is varied with P fixed. Of course, the multiple extrema show 
up here as well, as evidenced in Figure 12, along with some additional nontrivial de- 
pendence of the optimal on A. By comparing these bifurcation curves for different 
P, we see that the interval of multiple extrema moves to larger A as P increases (and 
vice versa). Furthermore, as A increases for fixed P, 0(P) at optimal can undergo 
abrupt increases, associated with the firing of an additional spike (Figure 13); the 
interesting structure of the BVP solution curves appears to be related to these events. 



5 LIF model with continuous input 

With the continuous input y (t) given in equation (28), we used variational methods to 
find values that yielded extremal values of 0(P) for the 0 model (27), given a final 
time P > 0. Such methods are not available for the LIF model due to the discontinuity 



Fig. 11 Time course of 6 for 
A = 8, P = 10, and various ^ 
values as labeled, which yield 
extremal values of ^(10). 
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Fig. 12 e{P) versus A for 
P = 10.5 (solid) and P = 13 
(dashed), along curves of 
extrema obtained by solving the 
BVP. The two curves coincide 
on the low, horizontal branch 
near 6{P) = 5. For each P, note 
the existence of multiple 
extrema for a small interval of A 
values as well as the occasional 
abrupt increases in 6, 
corresponding to spike addition. 




imposed by its reset condition. Nonetheless, we can perform some direct analysis of 
the dependence of firing in the LIF model on the parameter fi. In particular, although 
we will not specify an optimal fi, we will analytically establish some results about 
how spike times depend on ^, including the fact that the number of spikes saturates 
as p grows, which is consistent with [13] and contrasts with the non-monotone de- 
pendence of 0(P) on seen in the previous section. Similar methods can be applied 
to attain analytical results for the theta model, and these are briefly discussed in the 
Appendix; in particular, these show that spiking is lost when becomes sufficiently 
large. 

To perform this analysis for the LIF model, we make a fairly strong approximation. 
If a spike is fired at time Ta and the next spike occurs at time T^, then on the time 
interval (7^, 7),], we approximate y(t) hy the time average 



Tb - Ta J J 



Tb 



(34) 



Fig. 13 Time courses of 6 for 
optimal values of ^ found for 
P = 10.5 and various A values 
as labeled. Note that 6'(10.5) 
increases by slightly more than 
271 as A increases from 4.5, with 
6'(4.5)^ 5.052, to 10.5, with 
6'(10.5) ^ 11.771, and from 10.5 
to 16.5, with 6'(16.5) ^ 18.123, 
corresponding to the abrupt 
jumps in Figure 12. 
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with 

v^ = I-v-y{v-E), te (Ta, Ttl (35) 
Note that by definition, v(Tj~) = 0 and v(Tt) = 1, such that 

T T - C 

^ Jo I -v-y(v- E) 

(36) 

1 I^yE 

= In = . 

1+K I^yE-l-y 

If we assume that i;(0) =0 and fix a positive integer n, then we can try to solve 
equations (34) and (36) with Ta = 0 for a pair of positive numbers that we label 
(ri, Yi), such that the solution of equation (35) satisfies v(Ti) = 1. Similarly, we can 
set Ta = Ti and solve for (T2> Ti,y2), and inductively, given (Tj, yj), we can solve 
for (Tj-^i , K;+i) until we find T^+i > P for some n. At that point, we would propose 
that Ti, . . . ,Tn are our approximate spike times in the interval [0, P]. In practice, 
to constrain the solution set of this system of equations, we assume that T^ = P, 
since numerical explorations show that, for fixed A and P , if y0 is tuned to maximize 
the number of spikes generated on t e [0, P] by the LIF model with input given by 
equation (28), the final spike time is indeed generally close to P. Thus, for fixed A, 

P, instead of solving iteratively, we attain candidate approximate spike times by 
simultaneously solving n copies of equation (34) and n copies of equation (36) for 
unknowns [Ti, . . . , Tn-i, yi, • • • ,yn} with a final spike time T^ = P. 

We next study constraints on these approximate spike times {Ti, . . . ,Tn}. Specifi- 
cally, we have the following result. 

Proposition 5 For fixed A, P, there exists pi > 0 such that if p < pi, then there 
are no spikes (as defined above). If p > pi, then there exist 0 < < t(p) such 
that all spikes lie in [l(P),t(p)]. The function is monotone decreasing in p 
and is bounded above by There also exists P2 > such that the function t(p) 
is monotone increasing for p G (pi, ^2)^ achieves its maximum at p = p2, <^nd is 
monotone decreasing with I3t(l3) > 2 for p > 132. The value of ^2 is given by the 
minimal positive solution of fit (P) = 2. 

Proof First, note that spikes can only occur if Ap^te~^^ > -^^y . Clearly, the equation 

Aphe-^' = (37) 

has (Pi, ti) = ( ^(^"L^) ^ ^^eli-i) ^' ^^^^ = 1, as one solution. A brief calcula- 

tion verifies that equation (37) has no solutions for < fii, while for > fii, there 
are two solutions, (P,t_(P)) in {fit < 1} and (P, i(P)) in {fit > 1}. For fixed all 
spikes must lie in the time interval [l(fi), t(fi)]; note that, in particular, Ti is bounded 
below by l(fi), while our assumption that Tn ^ P would limit us to parameter choices 
such that t(fi) > P, although this inequality need not be satisfied for this proposition 
to hold. 
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Now, consider a solution t(P) of equation (37). Differentiation of (37) yields 

dt _ t(Pt -2) 

while equation (38) itself yields 



(38) 



= pt\P) + t{P) = (39) 

Equations (38), (39) imply that l(P), Pl(P) are monotone decreasing, since 

Pm < 1 (40) 

for all P for which it is defined. Initially, equation (38) also implies that t(P) is 
increasing, since the curve (P,t(P)) lies in 1 < < 2 for y0 near Pi. However, 
Pt(P) also increases, by equation (39), until t achieves a maximum at P2 such that 
PlKPl) = 2. For p > P2, i(P) decreases, by equation (38), but pt(P) > 2 for all 
P > P2, since the curve pt = 2 has a negative slope in the (P,t) plane and dt /dp = 0 
at pt = 2. □ 

An example of the solution curves to equation (37), illustrating the results of 
Proposition 5, appears in Figure 14. The curves of possible spike times Ti{P) lie 
between t_{P) and t{P) described above. Of course, if P is small, the spiking that 
can occur in time [0, P] will be limited. With large P, so that the full collection of 
predicted spike times is realized, direct simulations suggest that the number of spikes 
increases with P (for example. Figure 15). While we will not prove that result, we 
can establish some properties of spike times in the model in the limit of large P : 



Fig. 14 Plot of the solutions to 

equation (37) for A = 10, 

E = 1.2, I = 0.7. The curves 

shown are r = 2//3 (red), 

t = l/^ (blue), t = t{^) (dashed 

black), and t = t_{^) (large 

dotted black). Note that 

t_{^) < l/P and that t(P) rises 

up from r = l/yS to its maximum 

att = 2/^ and then decays, in 

agreement with Proposition 5. 
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Fig. 15 Number of spikes generated by the LIF model with A = 100, p = 10 and (a) E = 1.2, I = 0.7; 
(b) £^ = 2, / = 0.3. The blue curves show results of direct simulations of the LIF model with continuous 
input y(t) while the red lines are at A/\n(E/(E — 1)). 

Proposition 6 As P ^ oo, 

m ^ 0, 
m ^ 0, 

Pt{P) 00, and 
where n is the total number of spikes fired. 



lOU 



\ 100 




Proof The first Hmit is an immediate consequence of the bound (40). The second 
Hmit follows from equation (37). Specifically, some algebra yields 

A{E - _ _ 

(1 - I)efi ' 

with the left hand side clearly converging to 0 as ^0 ^ oo. Moreover, the inequality 



Afi^^t > Afi'^te 



E-1 
yields 

/AP(E-ir 



Pt>ln 



1 



so yS/"^ 00 as ^ 00. 

It remains to establish the fact that number of spikes converges as becomes 
sufficiently large, as stated in the Proposition and seen in Figure 15. First, note that 
the difference between successive spike times goes to 0 as y0 ^ oo, since ^ ^ 0 in 
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this limit. Thus, the left hand side of equation (36) tends to 0 as ^0 ^ oo. Maintenance 
of the equality in equation (36) therefore requires that 



y ^ 00 as ^ 00, 

which is consistent with the fact that the maximum of y(t) blows up with fi. 
Combining equations (34) and (36) yields the equality 



(41) 



A[pTae-^^- +e-P'^ - fine 



-^n _^-^^Tti _ _y_ 



In 



I + yE 



and equation (41) implies that the limit of the right hand side equals \n{E/{E — 1)). 
Hence, we define a constant Ci satisfying 



A[\-Cie-^' -^-^1] =ln 



E-\ 



(42) 



A solution of equations (34), (36), with To = 0, is Ti = Ci/fi, yi = (A- ACie~^^ - 
Ae~^^)fi/Ci. Next, define a constant C2 satisfying 



A[Cie-^' + e-^' - (Ci + C2)e-^'~^^ - e'^'-^^] = In 



Addition of equations (42), (43), after rearrangement, yields 



A[(Ci + C2)e-^'~^^ + e-^'-^^] = A - 21n 



E-1 



(43) 



E-1 



andri = Ci/^,72 = (Ci+C2)/)g, 

y2 = [Ci^-^i + e-^' - (Ci + C2)e-^'-^^ - e-^'-^^]Afi/C2 

is another solution of equations (34), (36). 

Repeating this process, we can find a series of solutions to these equations (see 
Figure 16), such that for the jth spike, Tj^i = Ci)/fi, Tj = (J^Ui Ci)/fi, 



.i=i 



Afi/Cj 



and 



= A-jln 



E-1 



(44) 
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Fig. 16 Times between spikes 
fired by the LIF model with 
A = 100, p = 10, ^ = 5, 
^ = 1.2, 7 = 0.7. The blue 
symbols show the results of 
direct simulations of the LIF 
model with continuous input 
y{t) while the red circles show 
solutions of equations (34), (36) 
obtained using the serial 
procedure described in the text. 
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We denote the total number of spikes by n(P) and estimate the final spike time 
r^(^) for any fixed pbyt(P), which gives Yll=i = Pt(P) ^ sls P ^ oc. Thus, 
from equation (44), it follows that lim^^oo ^ — ^(P) ln(E/(E — 1)) = 0, which yields 

n(p)^A/ln(^^^^ as ^^00, 

as desired (Figure 15). □ 



6 Discussion 

Summary and modeling issues 

We have considered how certain constrained, positive inputs should be timed to yield 
maximal numbers of spikes in the LIF and theta models for neurons. In both models, 
we have considered parameter regimes in which inputs must be above a threshold to 
elicit a spike. Thus, when each model is subjected to a train of discrete inputs of a 
fixed total magnitude, it is possible that maximal firing is attained by a critical kick 
strategy of giving just enough input to push the model trajectory above this threshold 
and then giving the minimal kick after each spike that achieves threshold clearance 
again. Aside from this possibility, we have analytically identified which other strate- 
gies could possibly be optimal for each model. Defining 8(g) as the magnitude of the 
change in input from the level g at the firing of one spike to the level g — 8(g) at 
the firing of the next spike, we found for the LIF model that our analytical approxi- 
mation to 8 (g) could be monotone decreasing in g or not, depending on the sign of 
E -\- 1 — 2EI An each case, we present an optimal strategy. If £^ + / — 2£^/ > 0, the 
optimum is a big kick strategy in which all available input is provided immediately. 
If E -\- 1 — 2EI < 0, if the minimum of 8(g) occurs at ^ = ^o, and if the amount of 
input available exceeds then the optimum consists of an initial kick of size ~ 
followed by kicks of size ^ 8 (go) after each spike, until the input is depleted. This di- 
chotomy of possible outcomes may be unexpected in light of standard intuition about 
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the LIF model, which emphasizes the power of big, synchronized input kicks to in- 
duce firing. Furthermore, the definition and analysis of 8(g) is itself novel, and its 
non-monotonicity in g means that there is a different effective dissipation of inputs at 
different input strengths, defined relative to the time it takes to progress from repolar- 
ization after a spike to firing of the next spike. This idea, that the power of a synaptic 
input depends on more than the individual input's amplitude, decay rate, and arrival 
phase, is often neglected in neuroscience studies and is an important observation of 
our work. 

Unlike the LIF model, our approximation to 8(g) always has a unique local mini- 
mum for the theta model, and we can compute its location directly. Assuming that this 
minimum corresponds to an initial condition that really does result in a spike, there 
is a clear optimal strategy of providing critical kicks to keep g near the minimum of 
5(g), as in the LIF case with E -\- 1 — 2EI <0. The existence of this minimum likely 
is related to the presence of an internal peak in the phase response curve (PRC) for 
the theta model, which past authors have suggested would represent an optimal phase 
for input timing [16], although PRC theory applies to oscillators receiving weak in- 
puts, while we consider strong inputs in the excitable regime. Heuristically, this min- 
imum may arise because of the shape of the ^-nullcline in the (0, g) plane. If g is on 
the small end of the spiking range, then the progress of 0 towards spike threshold is 
slowed by the proximity of the trajectory to the ^-nullcline, allowing g to drop signifi- 
cantly from one spike to the next. If this idea is correct, then we expect non-monotone 
8(g) to occur for models where voltage can be significantly slowed between spikes, 
without preventing the firing of the next spike, when g is on the lower end of the 
spiking range. In summary, for both the LIF and theta models with discrete inputs, 
the maximal number of spikes is attained for one of two strategies, either a critical 
kick strategy based on a minimal input threshold or a second big kick or critical kick 
strategy that we have identified from analysis of 8(g). The latter applies equally well 
in the case when the models are intrinsically oscillatory in the absence of inputs. 

We also show how to estimate analytically the number of spikes resulting from 
any input time course in both models, obtaining results that compare well with direct 
numerical simulations. Further, we establish an optimum value for the dependent 
variable of each model, v and 0 respectively, where each input should be delivered 
in a critical kick strategy to achieve maximal subsequent spike output. In both cases, 
this value is a stable critical point for the intrinsic dynamics of the model. 

When the input to each model is continuous, the optimization problem we consider 
is how to adjust the shape of the input, using a parameter that does not affect the 
input's overall magnitude, to achieve maximal spiking in a fixed time interval [0, P]. 
Derivation and solution of a BVP yields values of eliciting extremal numbers of 
spikes for the theta model, and numerical simulations show whether these extrema 
are maxima or minima. We find that the number of spikes fired increases and then 
decreases again as increases, corresponding to faster rise and decay of the input 
function and a larger maximal input. The existence of an interior local optimum for 
P is consistent with the non-monotonicity of 8(g) for the theta model in the discrete 
input case, with both observations pointing out that delivery of input in a stronger, 
faster way may reduce the resulting number of spikes for the theta model. Eventually, 
the number of spikes saturates, remaining invariant under additional increases in fi. 
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The LIF model with a continuous input shares this saturation property, as we prove 
in Proposition 6. For the LIF model, however, unlike the theta model, the number of 
spikes increases monotonically with and hence with the degree to which the input 
is concentrated in time, based on numerical simulations. This finding is consistent 
with previous analysis [12, 13] showing that synchronous input to the LIF model 
yields more spiking than inputs that are spread out. Indeed, in light of this set of re- 
sults, it is interesting that we do not always find this behavior for the LIF model with 
discrete input kicks, in the case where E -\- 1 — IE I < 0. This disparity in results for 
the LIF model points out that details of how synaptic inputs are modeled can influ- 
ence model dynamics in significant ways. The specific differences between optimal 
input patterns for the LIF and theta models that we highlight represent novel find- 
ings about the relationship between these models, while other differences have been 
demonstrated in earlier work [12]. Both of these models have Type I behavior [10, 
17], including the ability to fire spikes at arbitrarily low frequencies, yet clearly there 
are differences in their dynamic properties, including their responses to inputs. Such 
subtleties point out that classifications of models and neurons into gross categories, 
such as integrators and resonators, often need additional refinements to capture the 
diversity of neuron dynamics. 

Analytically, we found an interval of times during which all spikes must occur for 
the LIF model with continuous inputs, and we characterized the interesting depen- 
dence of the endpoints of this interval on ^. The results of applying similar methods 
to the theta model are also discussed in the Appendix. As ^ increases, all spike times 
converge to 0, yet the input becomes strong enough to elicit increasingly more spikes 
(up to some level). We were not able to exploit our approach to prove that the number 
of spikes increases monotonically with however. More specifically, we identified a 
minimal value Pi such that > fiiis necessary for spikes to occur. One idea for prov- 
ing monotonicity would be to seek a sequence of values fii < ^2 < ^3 < • • ' such that 
P > Pn^^ required for the firing n spikes (note that this ^2 differs from the ^2 used in 
Section 5). However, we were able to derive an equation for fii because we have an 
analytical formula for a minimal level of input needed for at least one spike to occur, 
and we do not have such an expression for subsequent spikes. Nonetheless, it is pos- 
sible that the band structure established in Section 2 could be exploited to abstractly 
establish results in this direction that might lead to a proof of monotonicity. 

The overall utility of the techniques presented in this paper depends in part on 
their generalizability. Note that the models we consider have relatively few parame- 
ters, and our results do not depend on particular parameter choices as long as they 
render the neuron excitable. For our numerical examples, we have chosen various 
parameter combinations to try to illustrate different parameter regimes. In some dis- 
crete input cases, we did choose values that appear to be too small to represent fast 
AMPA-mediated synaptic excitation, because certain differences in outcomes across 
strategies are most clearly evident with small . It is important to keep in mind, how- 
ever, that our models are dimensionless, that there are also slower NMDA-mediated 
excitatory synaptic currents, and that what we have represented as a slow synaptic 
decay could also arise from a long membrane time constant in a postsynaptic neuron 
or from the gradual arrival of many inputs from a presynaptic neuron population with 
a slowly down-ramping firing rate. Aside from parameter variations, our methods for 
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estimating numbers of spikes are also quite general across different discrete input 
patterns, for the models we consider. For these models, our methods could likely be 
adapted to other optimization problems, such as tuning inputs to achieve regularity of 
interspike intervals, spiking within some range of rates, or spiking at particular times 
[7] or finding minimal inputs to generate particular spike patterns [8]. In each of these 
problems, we would obtain predictions about which features of the timing and size 
of inputs within a discrete input train yield which spiking patterns. 

We have assumed that the neurons receiving inputs are silent in the absent of in- 
puts. As noted above, our techniques for the discrete input case would generalize to 
models with nonzero intrinsic firing rates, and in particular the band structure par- 
titioning phase space into initial conditions corresponding to different numbers of 
spikes fired before inputs are depleted would carry over analogously. Indeed, the ex- 
istence of a minimal level of input needed for spiking becomes irrelevant during time 
periods in our analysis when the input is well above that level, such as during the 
application of critical kicks based on a local minimum of 8(g). The presence of in- 
trinsic oscillations complicates the case of continuous inputs, because the impact of 
an input will depend on the phase at which it starts, as in recent phase resetting anal- 
ysis for bursting models [18, 19]. Furthermore, our results about constraints on input 
times during which spiking can occur in the LIF model with continuous inputs clearly 
depend on the absence of intrinsic firing. 

Unlike the models that we have considered, many other neuronal models have 
two- or higher-dimensional intrinsic dynamics. Reductions into fast and slow subsys- 
tems (for example, [20, 21]) or other reductions [22] offer the possibility of extract- 
ing subsystems from these models on which a similar analysis to ours, including the 
partitioning of an appropriate phase space into bands associated with certain spike 
numbers in the discrete case, could be performed. It is also possible that results could 
be obtained after reduction to a firing time map (for example, [23]). Moreover, ana- 
lytical techniques can be used to reduce general oscillator models of arbitrary dimen- 
sion, subject to weak inputs of a prescribed time course, down to forced scalar phase 
equations [24] . Optimization techniques could be applied to tailor the forcing terms 
in such equations, within certain constraints, but of course this option is only avail- 
able if the intrinsic neuronal dynamics is oscillatory and the inputs are weak. Direct 
simulations can also be done to begin to examine how 8(g) varies with g for par- 
ticular higher-dimension models. Some preliminary simulations yielded a monotone 
8(g) for a few conductance-based models (for example, [25-27]) in certain param- 
eter regimes, as well as a non-monotone 8(g) for a particular scaling of a reduced 
Hodgkin-Huxley model ([27], Type B~), but a thorough exploration of such models 
remains for future work. 

Neuronal relevance and related issues 

Given that we can identify an optimal input structure, it is important to ask whether 
the similarity of actual input streams to the proposed optimum can be checked experi- 
mentally and whether such an optimal input could actually be realized by a network of 
neurons. As for the former question, membrane potential dynamics can be recorded, 
including identification of changes in potential associated with synaptic inputs, and 
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associated synaptic input features can be estimated [28], so it appears that the relevant 
experimental techniques are indeed available. 

As for the latter, at least certain features of the optimal inputs we have identified do 
appear to be achievable. In the discrete input case, consider first a pool of presynaptic 
neurons providing inputs with the same decay rates, with some distribution of input 
amplitudes. This pool could be large or small, depending on the brain area involved. 
The overall input to a postsynaptic cell could be tailored by adjusting the relative fir- 
ing times of these presynaptic cells. For example, a critical kicks strategy not based 
on a minimal input threshold could be achieved if a large initial input was given, fol- 
lowed by a sequence of appropriately timed smaller inputs. This input pattern could 
be realized either through a synchronized input from a group of neurons followed 
by later inputs from a subset of these neurons or from a distinct group of neurons, 
perhaps made active by the initial large input as well. Short-term synaptic depression 
could play a role in this tapering of inputs within the input train. The presence of 
short-term plasticity would alter the set of input patterns that could be provided, but 
for a fixed form of postsynaptic dynamics, only certain input patterns would be near- 
optimal, and if short-term plasticity promoted these, then the loss of the capability 
to produce non-optimal patterns would be irrelevant. Similarly, neuromodulators also 
can modulate synaptic dynamics, in a way that differs across different cell types [29- 
31]. Although we have neglected short-term plasticity and neuromodulation here, 
considering their effects on optimal input strategies could be an interesting direction 
for future work. Moreover, if a given presynaptic neuron repeatedly fires in a certain 
timing relationship with a postsynaptic neuron under an optimal input pattern, then 
spike timing dependent plasticity (STDP) could also become relevant. Generally, the 
properties of the inputs received by a given postsynaptic neuron, or the membrane 
potential variations in that neuron that are induced by inputs, can vary across differ- 
ent situations in real neuronal networks; indeed, these inputs need not come from the 
same presynaptic source in different states [32]. This observation is consistent with 
the idea that input patterns could be tailored to achieve different functions. 

Given these arguments in favor of the idea that inputs to neurons can indeed be var- 
ied in biologically reasonable ways, our work leads to several biological conclusions 
and predictions. First, the pattern of inputs needed to maximize firing will depend 
on a neuron's intrinsic dynamics. This point has been discussed in previous theo- 
retical work in the weak or noisy input limit (for example, [27, 33, 34]) as well as 
in some past experimental work [35, 36], and we expect this principle to hold quite 
generally. Second, there may exist an optimal timing following a spike at which the 
input strength needed to elicit an additional spike is minimized. Again, although we 
consider excitable neurons receiving strong inputs, this idea is related to past work 
on PRCs for oscillators receiving weak inputs; this idea differs from the concept of 
resonance, however, in that the optimal timing would not be solely determined by a 
postsynaptic neuronal intrinsic frequency. Related to this idea of optimal timing, we 
might predict that neurons' afterhyperpolarization time courses would differ across 
brain areas that receive input streams with different characteristics (cf. [37]); that 
is, the intrinsic dynamics of afterhyperpolarization and network input characteristics 
could have co-evolved to achieve some form of optimality. Third, we would specu- 
late that perhaps background input levels could be tailored to keep neurons in active 



^ springer 



Page 32 of 35 



Wang et al. 



brain regions in a regime where 8(g) is small and hence the barrier to spiking is low- 
est; attention could even be related to the selection of such a state. Past results have 
shown that input fluctuations can establish a regime of high spike-time reliability 
[38, 39] or high sensitivity of firing frequency to input current strength [37, 40, 41], 
and tuning background input levels to facilitate or suppress postsynaptic activity is at 
least biologically plausible. Finally, we would also predict that excitatory postsynap- 
tic potential time courses themselves would differ in different neuron types, activity 
states, developmental stages, and brain areas, as has been seen in experimental work 
[30, 42], in a way that is interrelated with the dynamic properties of the postsynaptic 
neurons involved. 

Admittedly, the precise optimal input patterns that we have found in the discrete 
case do depend on certain aspects of the postsynaptic neuronal dynamics that might 
render them hard to achieve. For example, for the LIF model, we found that optimal 
inputs arise when v = I, but it is not clear how this information would be available to 
the input source. Clearly, we are considering a situation where the input source is not 
directly impacted by the firing of the postsynaptic neuron receiving the inputs; recur- 
rently connected networks are outside the scope of this work. Moreover, stochastic 
effects could alter the details of any target spike train and synaptic release pattern. 
It is reasonable to expect that stochastically perturbed versions of the optimal input 
patterns identified for deterministic models would provide the optimal response dis- 
tribution in the presence of weak noise, but careful investigation of this issue remains 
for future work. 



Similar methods to those applied to the LIF model with continuous input can be 
used for the theta model with continuous input. Here we concisely summarize the 
approach and the main results. Again, we replace y(t), given in equation (28), with 
its time average, and thus we consider the equation 



on the time interval (Ta,Tb) from one spike to the next. As 0(Ta) = — tt and 
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= 1 - cos6> + (/? + + cos^) 



with 




n-Ta 




= n-Ta. 
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Combining these equations yields 



TC 



,2 



(45) 



-b(Tt-Ta), 



a 



and spike times Ti and associated values yi are found by solving these equations 
repeatedly. 

Lower and upper bounds t(fi) on firing times are estimated from the solu- 
tions of 



since y (t) > —b > 0 is required for spiking to start. Technically, a single final spike 
can be fired if y rises above —b and then falls below it, however, so the upper bound 
here is not precise. 

As with the LIF model, constraints on spike times and changes in spike times with 
can be explored with these ideas. As one example, consider time Ti(fi) when the 

first spike is fired. It can be shown that t(fi) ^ 0 as ^0 ^ oo, and hence Ti(fi) 0 

as y0 ^ 00 as well, if it exists. This result yields 



but for Ta = To, Th = Ti, the left hand side of equation (45) becomes A — 
ApTie~^^^ — Ae~^^^ < 2A, a contradiction. Hence, no first spike can occur, and 
we conclude that spiking in the theta model with input (28) is lost as ^ oo. 
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